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SUMMARY 


New, improved curve fits for the transport properties of equilibrium 
air have been developed. The curve fits are for viscosity and Prandtl 
number as functions of temperature and density, and viscosity and thermal 
conductivity as functions of internal energy and density. The curve fits 
were constructed using Grabau-type transition functions t model the 
transport properties of Peng and Pindroh. The resulting curve fits are 
sufficiently accurate and self-contained so that they can be readily 
incorporated into new or existing computational fluid dynamic codes. The 
ranges of validity of the new curve fits are temperatures up to 15,000 K 
and Amagat densities (p/po) from 10~5 to 10*. 
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SYMBOLS 


e specific internal energy, m z /s z 

h specific enthalpy, m z /s 2 

k coefficient of thermal conductivity, J/K-m-s 

p pressure, N/m 2 

Pr Prandtl number 

T temperature, K 

p dynamic viscosity, kg/m-s 

p density, kg/m 3 

Subscript: 

o reference conditions of 273.15 K and density of 1.243 kg/m 3 
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INTRODUCTION 


The realm of high-speed flight requires the knowledge of the 
transport properties of equilibrium air. Tie transport p operties 
for air In equilibrium can be found in the literature in the form 
of tables (refs. 1 and 2). However, in performing computational 
fluid dynamics calculations, the resulting table-lookup process is 
often time consuming and cumbersome. As a consequence, s 1 e 
closed-form equations that can be readily incorporated Into 
existing numerical codes are desirable for the various transport 
properties . 

Vigneron (ref. 3) developed simplified approximations for p = 
y(e,p) and k = k(e,p), for the above reasons. The data for these 
curve fits were derived from Hansen (ref. 1). One of the major 
shortcomings of Vigneron's curve fits is the relatively la-ge 
error, especially in the curve fit for k = k(e,p). The primary 
objective of the present research effort was to develop a curate 
curve fits for the following transport properties: 

V = P(T>p) 

Pr = Pr(T,p) 

V = P(e,p) 

k = k(e,p) 

The curve fits are based on the transport properties computed 
by Peng ard Pindroh (ref. 2) and have the same ranges of validity, 
namely, temperatures up to 15,000 K and Amagat densities (p/p 0 ) 
from 10~5 to 10 1 . 
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SOURCE OF EQUILIBRIUM AIR TRANSPORT PROPERTIES 


It is appropriate at this point to briefly describe the approach used 
by Peng and Pindroh to calculate the transport properties. They chose a 
nine species (N 2 , O 2 , NO, 0, N, 0 + , N + , N0 + and e") model )f air 
as an adequate representation for the transport properties under 
consideration. The non-uniformities of the mass-velocity temperature, 
and composition of the gaseous system result in the trans )rt processes of 
momentum, energy and mass transfers. The interaction for is between the 
gas molecules determine the transport mechanism. 

The relationship between the transport properties ano the 
intermolecular forces is given by the Chapman-Enskog method. In this 
method one solves for the transport coefficients In terms of a set of 
collision integrals. It is also possible to express the transport 
coefficients in terms of rigid sphere cross-sections. This was the 
procedure adop- ed by Peng and Pindroh, who used experimentally determined 
cross-sections. The Chapman-Enskog method is valid for situations when 
only binary collisions are important, and where the gradients of the 
physical quantities are small. 

The Prandtl number computed by Peng and Pindroh and correlated in 
this report is based on the frozen property concept, i.e., it does not 
include tie reactive component due to diffusion as discussed by Hansen 
(ref. 1). 


CONSTRUCTION OF THE CURVE FITS 

Figure 1 shows the variation of p/p 0 with temperature at various 
constant tensities. These curves exhibit segments of linear or quadratic 
functions successively connected by transition functions which are 
asymptotic at both ends, and include points of inflection Reference 
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4 outlines in detail how these curves can be systematically approximated 
by exponential transition functions using the technique of Grabau 
(ref. 5). Following the method of Grabau, one has a choice of two 
kernel transition functions. The first is the Fermi-Dirac function 

, L_. (l) 

1 + exp(mx) 

which represents a transition with an inflection between the levels 
0 and 1. This transition function will be referred to as an odd 
transition function. The second type of transition function is the 
kernel of the Bose-Einstein distribution function 
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( 2 ) 


1 - exp(mx) 

which provides a transition without a point of inflection, and is 
known as an even transition function. In both these transition functions, 
the sign and magnitude of the constant m determine the direction 
and rate of the transitions, respectively. 

These kernel transition functions can be modified to nodel transitions 
in two independent variables. The kernel of an odd trans tion function 
in three dimensions is 

1 (3) 

1 + exp(a Q + ajX + a 2 y + a 3 xy) 

where x and y are the two independent variables. Equation (3) is 
essentially an alternate form of 

1 (4) 

1 + expLm(x-x o )(y-y o )] 

However, equation (3) is more convenient for determining the values 
of the constants a^ and a^ as dictated ly the behavior imposed on 
the transition function. The general technique of determining the 
values of these constants is as follows. The boundaries of the transition 
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in the directions of the two independent variables are x a - x i x b 
and y c s y sy d . If f^x.y) and f 2 (x,y) are the two surfaces limiting 
the transition function f(x,y), then 

f?(x,y) - f i (x,y) 

f(x,y) = f,(x,y) + ; — (5) 

1 ± exp(a Q + ajX + a ? y + a 3 xy) 

In order to ensure an accurate and smooth transition from f^(x,y) 

to f 2 (x,y), it is necessary for the quadratic expression (a Q + a^x 

+ a 2 y + a 3 xy) to behave as follows. At the lower left corner point 

(x , y ) the quadratic expression should have a large positive value 
s c 

so that f (x,y) ^f^x.y). At the upper right orner point (x b> y d ), 
the quad’atic expression should have a large negative value in order 
to ensure that f(x,y) - f 2 (x,y). At the midpoints of the left and 
right boundaries, Cx a ,(y c + y d )/2] and [x b ,(y c + y d )/2], respectively, 
the quadratic expression should be zero so that 

f i (x,y) + f ? (x,y) 

f(x,y) =— -i - 

2 

These conditions yield the following four linear equations: 


+ a l x a + a 2 y c + a 3 x a y c = +m 

(6) 

+ a lXh + a 2 y d + a 3 x b y d = -m 

(7) 

+ a l x a + a 2 (y c + y d )/2 + a 3 x a (y c + y d )/2 = 0 

(8) 

+ a i x b + a 2 (y c + y d )/2 + a 3 x b (y c + y d )/2 = 0 

(9) 


where m is a positve constant (typically, 10 sm <20) chosen such 
that exp (m) and e>p (-m) do not yield overflow or underflow conditions, 
respectively, on a computer. The constants a Q through a 3 can now 


7 


be obtained in a straightforward manier from the system of four linear 
equations in four unknowns (eqs. (6) - (9)). 

The above method of obtaining the Grabau-type transition functions 
proved quite accurate in ensuring a negligible mismatch in the dependent 
variable over the boundaries of adjoining sub-regions. Graphical 
inspection of each of the transport property surfaces al ows one 
to identify the various sub-regions a: well as the exact locations 
of the necessary Grabau transition furctions. It is a merit of this 
stepwise method of constructing empirical equations that any part 
can be removed for corrections without disturbing the sur ace approximation 
as a whole. 


EQUATIONS OF THE CURVE FITS 


The curve fits for the various transport properties were constructed 
using Grabau-type transition functions, as described previously. 

The general form of these curve fits can be written as 


, x , , " M*’ 7 ) 

z x,y f^x.y) + j ± exp ( m ^ + + + 


( 10 ) 


whare, in general, 

^(x.y) = pj p 2 x + p 3 y + p 4 xy + p 5 x 2 + p 0 y 2 + p ? x 2 y + p g xy 2 

+ P<) x3 + Pi 0 y (ID 

and 

f 2 ( x »y) = Pn + p 12 x + p^y + Pi 4 xy + p 15 x2 + p 16 y 2 + P 17 x 2 y 

+ Pl8 xy2 + Pl9 x3 + P20 y3 (12) 
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The positive sign in the denominator corresponds to an o< d transition 
function while a negative sign corresponds to an even trnsition 
function. The coefficients m Q through m^ in the denominitor of the 
transition function in equation (10) were determined by the technique 
outlined in the preceding section. The coefficients Pj through PgQ, 
in equations ill) and (12), were determined by the actual curve fitting 
of the dita from Peng and Pindroh (ref. 2). The exact location 
and number of these data points over the curve fit domain determined 
the accuracy of the curve fits. The points were clustered near the 
boundaries of the domain and the mid-region of the transition in 
order to ensure continuity across the boundaries and accuracy within 
the doma n. The data of Peng and Pindroh were fitted to he equations 
of the curve fits by the method of least squares. The FORTRAN subroutine 
employed to perform the least squares curve fitting of the data uses 
a multiple linear regression technique (ref. 6). 

Pen} and Pindroh tabulate the three transport properties p, 

Pr and k in increments of 100 K for temperatures from 500 K to 15,000 K 
along 13 constant density lines. In order to generate data at intermediate 
points for the least squares curve fitting, an interpolation procedure 
was necessary. Two IMSL subroutines IBCCCU and IBCEVL (ref. 7) were 
used to perform bicubic spline interpolation of the data. The interpolated 
data v/as found to be smooth and was therefore, used with good confidence. 
The general form of the curve fit for each transport property is 
described below. 


9 


1 LI M(T,p) 

The general form of the equation used for the correlation 
H = u(T, o) was 


H/W 0 


il l + a 2 X + a 3 Y + a 4 XY + a 5 x2 + a 6 Y2 + a 7 X 2 Y + a g XY 2 + a 9 X 3 

+ a 10 Y3 + < a ll + a 12 X + a l 3 Y + a l4 XV + a 15 )|! + a 16 V! * a j7 
+ a 18 XV a + a 19 x 3 + a ?0 ¥ 3 )/[l t exp(a 2 , + a 22 X t a 23 Y + a 2 „ 


X 2 Y 

XY)] 

(13) 


where X = T/1000 and Y = log 10 (p/p o ). The units for T are K and 
the units for p are kg/m 3 . It should be noted that not all the terms 
appearing in equation (13) are used over the complete range of T 
and p. The units of p are kg/m-s. 

Pr = Pr(T, 0 ) 

The general form of the equation used for the correlation 
Pr = Pr(T,p) was 

Pr = bj + b 2 X + b 3 Y + b 4 XY + b g X 2 + b g Y 2 + b ? X 2 Y + b g XY 2 + b g X 3 

+ b 10 Y3 + (b ll + b 12 X + b 13 Y + b 14 XY + b 15 XZ + b 16 Y2 + b 17 X2Y 
+ b lg XY 2 + b lg X 3 + b 2Q Y 3 )/[l + exp(b 21 + b 22 X + b 23 Y + b£ 4 XY)] 

(14) 

where X = T/1000 and Y = logio(p/p 0 ). The units for T are K and 
the units for p are kg/m 3 . 
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it = M( e »P? 

The equation for the correlation p = p(e,p) was 


p/p Q = c x + c 2 Z + c 3 Y + c 4 YZ + c 5 Z 2 + c g Y 2 + c ? Z 2 Y + c g ZY 2 + c g Z 3 

+ c 1Q Y 3 + (c n + c 12 Z + c 13 Y + c 14 YZ + c 15 Z 2 + c 16 Y 2 + Cj 7 Y 2 Z 
+ CigYZ 2 + c 19 Z 3 + c 2Q Z 3 )/[1 ± exp(c 21 + c 22 Z + c 23 Y + c 24 YZ)] 

(15) 


where Z = logig(e/e 0 ) and Y = logio( p / f b)- The units of u are kg/m-s, 
the units for e are m 2 /s 2 , and the units for pare kg/m 3 . 

k = k(e, p) 

For the correlation k= k(e,p), the general form of the equation 
used was 


k/k Q = d 1 + d 2 Z + d 3 Y + d 4 YZ + d g Z 2 + d g Y 2 + d ? Z 2 Y + d g ZY 2 + d g Z 3 

+ d 10 Y3 + (d ll + d 12 Z + d 13 Y + d 14 YZ + d 15 Z2 + d 16 Y2 + d 17 Z2Y 
+ d lg ZY 2 + d ig Z 3 + d 2Q Y 3 )/[l ± exp(d 21 + d 22 Z + d 23 Y + d^YZ)] 

(16) 

where Z = logio( e / e o) and Y = logio (c/po)- ^he units for e are 
m 2 /s 2 , the units for pare kg/m 3 , and the units of k are Joule/ K-m-s. 


RESULTS AND CONCLUSIONS 


New, simplified curve fits for the transport properties of equilibrium 
air were constructed using the methodology described in the f recedi lg 
sections. Comparisons of the curve fits p = p(T,p), Pt = Pr(T,p), 
p = p(e,p) and k = k(e,p) with the data of Peng and Pindroh (ref. 

2) are shown in figures 2 to 5. The following procedure was employed 
in generating the curve fits for p = p(e, p) and k = k(e,p). For 
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°n/i 
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16.5 
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0.500 0.875 1.250 1.625 2.000 2.375 2.750 3.12S 3.500 

L°gio(e/ e 0 ) 

Figure 4. Comparison of curve fits for p = p(e,p) 




a given pressure p and density p, the NASA R6AS program (ref. 8) 

(which gives the thermodynamic properties of equilibrium air) was 
used to obtain the corresponding h and T. The given p and computed 
T were then used to obtain p and k from a bicubic spline interpolation 
of the data of Peng and Pindroh. The corresponding value of e is 
given by e = h - o/p . This same procedure was used for the present 
comparisons. 

The comparisons are presented graphically :o provide a qualitative 
overview of the accuracy of the curve fits. Hovever, as figures 
2 to 5 indicate, these graphical comparisons are restricted to points 
lying on 7 constant-density lines ranging from 10 -5 to 10 1 Amagats. 

In order to ensure the validity and accuracy of the curve fits across 
the entire domain, a more comprehensive accuracy test was carried 
out. The curve fits were compared with the data of Peng and Pindroh 
for relative accuracies at approximately 20,000 data points. These 
test points were chosen to span the entire density range from 10" 5 
to 10 1 Amagats and temperatures varying from 500 K to 15,000 K. 

The results of these comprehensive accuracy checks are presented 
in tables 1 to 4. The first column in the tables represents the 
percentage e* ror in the comparison of a property generated by bicubic 
spline interpolation of the data of Peng and Pindroh and a curve 
fit. The second co umn contains the ,>ercen:age of test points which 
are in e ’ror by an amount greater than that indicated in column 1. 

The accuracies of the curve fits for p = p(T,p) and Pr = Pr(T,p) 
are good to within 4 percent of the data of Peng and Pindroh. The 
somewhat higher percentage errors in the two curve fits with e and p 
as independent variables can be attributed to the steep gradients 
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Table 1. Accuracy of u = h(T,p) 


Error, 

percent 

Percentage of data points 
with higher error 

0.5 

26.19 

1.0 

9.01 

2.0 

1.05 

3.0 

0.14 

4.0 

0 

5.0 

0 

6.0 

0 

7.0 

0 

8.0 

0 

9.0 

0 

o 

o 

r-H 

All 

0 


Total number of data points = 35,211 

Maximum error = 3.71 percent at T - 15,000 K and logjo( p / p o) = "5.0 
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Table 2. Accuracy of Pr = Pr(T,p) 


Error, 

percent 

Percentage of data points 
with higher error 

0.5 

13.10 

1.0 

0.77 

2.0 

0.05 

3.0 

0 

4.0 

0 

5.0 

0 

6.0 

0 

7.0 

0 

8.0 

0 

9.0 

0 

o 

o 

i-H 

All 

0 

Total number of data points = 

35,211 


Maximum error = 2.47 percent at T = 9350 K and 1ogio( p/pq) = 10~ 5 
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TABLE 3. Accuracy of m = p(e,p) 


Error, Percentage of data points 

percent with higher error 


0.5 

16.11 

1.0 

4.55 

2.0 

1.48 

3.0 

0.83 

4.0 

0.48 

5.0 

0.20 

6.0 

0.06 

7.0 

0.02 

8.0 

0.01 

9.0 

0.01 

o 

o 
«— 1 

All 

0 


Total nunber of data points = 18,794 

Maximum error = 9.22 percent at logio(e/e 0 ) = 3.28, T = 13,920 K and 

logio(p/p 0 ) = -5 -° 
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TABLE 4. Accuracy of k = k(e,p) 


Error, 

percent 

Percentage of data points 
with higher error 

0.5 

62. 6f 

1.0 

44. 72 

2.0 

28.40 

3.0 

17.40 

4.0 

4.14 

5.0 

2.58 

6.0 

1.69 

7.0 

1.04 

8.0 

0.62 

9.0 

0.37 

^10.0 

0.25 


Total number of data points = 18,794 

Maximum error = 15.27 percent at logio(e/e 0 ) = 2.65, T = 6,110 K and 
logio(p/P 0 ) = -5.0 
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in the surfaces which are a functional representation of the 
transport properties. The present approach of using Grabau transi- 
tion functions in conjunction with a least squares curve-fitting 
procedure cannot model the steep gradients to any greater accuracy. 
The curve fits for y = y ( e , p ) and k = k(e,p) were not compared with 
those developed by Vigneron (ref. 3) because the latter are based on 
the data of Hansen (ref. 1) which are substantially different from 
that of Peng and Pindroh. 

In conclusion, the new, simplified curve fits provide a fast and 
accurate means of obtaininj the transport properties of equilibrium 
air. They can be Incorporated into computational fluid dynamics 
computer codes in a straightforward manner. Their accuracies are 
sufficient to justify their use in both steady-state and time- 
dependent flow computations. 
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onnnoooooo 


APPENDIX — FORTRAN SUBROUTINE LISTINGS 


SUBROUTINE UGAS1 (T, RHO, MU) 

INPUTS FOR SUBROUTINE : 

T = TEMPERATURE, IN KELVIN 
RHO = DENSITY, IN KG/M**3 

OUTPUT : 

MU = DYNAMIC VISCOSITY, IN KG/M-S 

REAL T , RHO , MU , RHOO , Z , X , Y , GAS1 , GAS2 , JAS3 , GAS4 , 

*GAS5 ,GAS6 ,GAS7 ,GAS8 , GAS 9 
DATA RHOO/1 . 243/ 

X=T/1000 . OEOO 
Y= ALOG 1 0 ( RHO/RHO 0 ) 

IF ( (Y.LT.-5.EOO) .OR. (Y.GT.1.E00) .OR. (X.GT. 1 . 5E01 ) ) 
*WRITE (6,1000) RHO,T 
IF (T.GT.3O0.E00) GO TO 5 
MU=1.462E-O6*SQRT(T)/(1.0E0O+112.OE0O/T) 

RETURN 

5 IF (T.GT. 6O0E0O ) GO TO 10 

GAS1=4.13906E-01+2.16606E00*X 
GAS2=(1.30718E-05+7.44367E-05*X)*Y 
GAS3=(-5.45043E-02-1 .74550E-04*Y-1.15324E-01*X)*X*X 
GAS4=(2.43199E-05-2. 14485E- 05*X+3 . 03976E- 06*Y)*Y*Y 
Z=GAS1+GAS2+GAS3+GAS4 
GO TO 100 

10 IF(T.GT. 5500. OEOO) GO TO 20 

GAS1=4 . 8653102 E- 01+2. 1053953E+00*X 

GAS 2= ( -4. 4502862E-02+5. 062232 5E~02*X)*Y 

GAS3=(-2. 3327267E-01-1.019074E-02*Y+1.9685295E-02*X)*X*X 

GAS4= (2. 756468 0E -04+2. 7222564E~03*X+8 . 6649903E~04*Y) *Y*Y 

Z=GAS1+GAS2+GAS3+GAS4 

GO TO 100 

20 IF (T.GT. 1 . 05E04 ) GO TO 40 
IF (Y.GT.-2.5E00) GO TO 30 
GAS1=5.993881E01-1.698837E01*X 
GAS 2= (2.113989E01~5.130287EO0*X)*Y 

GAS 3= (1.7423 64 E00+2. 778705E-01*Y-5 . 201473E-02*X) *X*X 
GAS 4= (1. 63767 5E00-2. 460623E~01*X+1 . 671319E-02*Y) *Y*Y 
GAS 5= 4 . 43870 6E02-5.640044E01*X 
GAS 6= (7. 553438 E0 1~3. I77507E00*X)*Y 
GAS7=(2.078559EOO-7.210869E-02*Y-1.669889E-02*X)*X*X 
GAS 8= (5.507033E00+1. 196004E-02*X+2 . 199922E-01*Y) *Y*Y 
GAS9=EXP( 107 . 00E00-7 . 40E00*X+11 . 50E00*Y-0 . 41*X*Y) 

GO TO 90 

30 GAS1=3 . 5331 6E00+4. 9342 5E~01*X 

GAS2=(-4.06143E-01+2.10671E-01*X)*Y 

GAS 3= (7. 3493 4E -02-2. 6395 2E-02*Y-1. 26658E-03*X) *X*X 


ORIGINAL’ PAGE IS 

°R POOR QUALITY/ 
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GAS4=(-1 .30624E-01+2.74790E-02*X-2.79567E-03*Y)*Y*Y 

GAS5=2.00538E01-6.67992EOO*X 

GAS6=(1.03098E01-2.48845E0O*X)*Y 

GAS7=(7. 57575E-01+1 ,52736E-01*Y-2.91i22E-02*X)*X*X 
GAS8=(1.40864EOO-1.80169E-01*X+4.70062E-02*Y)*Y*Y 
GAS9=EXP(63.75EO0-7.976EO0*X+5.357E-Ol*Y+8 .333E-01*X*Y) 
GO TO 90 

40 IF (T.GT.13.0E03) GO TO 60 
IF (Y.GT.-4.0OEO0) GO TO 50 
GAS1=3 . 2488 5E02-5. 4635 9E01*X 
GAS2 = ( 2.86103E01 + 1.1196", E00*X)*Y 
GAS3=(4.1014lE00-2. 3842^ E-02*Y-1.08784E-01*X)*X*X 
GAS4=(4.44362E00+2.2490:E-01*X+4.59308E-01*Y)*Y*Y 
GAS'i = -4.50893E02 + 3.61004E01*X 
GASi.= (-1.46489E02 + 1.39297E01*X)*Y 
GAS/=(3.63567E-01-1.33686E-01*Y-2.80207E-02*X)*X*X 
GAS8=(-8.36798E00+8.40047E-01*X+1.77782E-01*Y)*Y*Y 
GAS9=EXP(2.447E02-1.874E01*X+4.856E01*Y-3. 723E00*X*Y) 
Z=GAS1+GAS2+GAS3+GAS4+ (GAS5+GAS6+GAS7+GAS8 )/( 1 . 0-GAS9 ) 
GO TO 100 

50 IF ( Y.GT. -2 . 5EO0 ) GO TO 80 

GASl=4.74364EOl-2.52946E00*X 
GAS2=(-3.40953E01+4.33761EOO*X)*Y 
GAS3=(4.19920E-02-1 .12842E-01*Y+9.09826E-05*X)*X*X 
GAS4=(-6.95452EOO+3.95061E-01*X-3.33072E-01*Y)*Y*Y 
GAS 5 = - 3 .457 58E02 + 6. 5481 2E01*X 
GAS6=(-3.77086E01+4.97501EOO*X)*Y 
GAS7=(-4.17677EOO-1.59916E-01*Y+9.0358E-02*X)*X*X 
GAS8=(2.01908EOO-8.41274E-02*X+2.61401E-01*Y)*Y*Y 
GAS9=EXF (-197. OEOO+14.6EOO*X-41.2EOO*Y+2.85EOO*X*Y) 

GO TO 90 

60 IF (Y.GT.-4.00E00) GO TO 70 
GAS1=4.53184E02-5.27482E01*X 
GAS2=(1.29609E02-9.91921E00*X)*Y 

GAS3=(2.07 504EOO + 1.90185E-01*Y-2.761C 9E-02*X)*X*X 
GAS4= (1. 2675 5E01-4. 8895 5E~01*X+ 4. 0942 8E-01*Y ) *Y*Y 
GAS5=-1.15162E02-4.02569E00*X 
GAS6=(-8.84578E01-2.17376E01*X)*Y 
GAS7=(-3.8511EOO-9.30247E-01*Y-2.00626E-02*X)*X*X 
GAS8= (-4.78794E01-4. 61508EOO*X-7 . 34409E00*Y ) *Y*Y 
GAS 9= EXP (76.82-2.29*X+15.08*Y-O . 4475*X*Y ) 
Z=GAS1+GAS2+GAS3+GAS4+ (GAS5+GAS6+GAS7+GAS8 )/ ( 1 . 0~GAS9 ) 
GO TO 100 

70 IF (Y.GT.-2.50E00) GO TO 80 
GAS.'. = 4 . 5228 9E 02-5. 5093 2E01*X 
GAS 2= (1.18987E02-9.05427EOO*X)*Y 
GAS3=(2.34316E0O+1.857E-01*Y-3.4O138E-02*X)*X*X 
GAS4=(1.13829E01-3.91996E-01*X+4.075E-01*Y)*Y*Y 
GAS5=-7.93738E02+1.27668E02*X 
GAS6=(-1.76755E02+2.93067E01*X)*Y 
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GAS7=(-6.18382E00-7.86123E-01*Y+9.67472E-02*X)*X*X 
GAS8=( 2. 67668EOl+l.O8936EOO*X+7.2087UEOO*Y)*Y*Y 
GAS 9= EXP (~63.33E00+3.33E00*X-16.67E00*Y+0.667E00*X*Y) 

GO TO 90 

80 GAS1=5 .05519E02-6. 60139E01*X 

GAS2=(1.20621E02-9.43025E00*X)*Y 
GAS3=(3.07622E00 + 1.93366E-01*Y-5 . 10125E-02*X) *X*X 
GAS4=(1.10339E01-4.07027E-01*X+3 . 59506E-01 *Y) *Y*Y 
GAS5=-5 . 1450 5E 02+6. 9101 6E01*X 
GAS 6 = ( - 1.15237E02 + 8.06602E00*X)*Y 
GAS7 = (~3. 12671 E0 0-1. 1347 3E-01*Y+4. 8532E~02*X) *X*X 
GASS = ( -8 . 7 54 53E00 + 1 . 69286E-01*X-2 . 57493E~01*Y) *Y*Y 
GASS =EXP( -156. 1EOO+9. 58E00*X~32. 3E00*Y+1 ,64E00*X*Y) 

90 Z=GAS1+GAS2+GAS3+GAS4+ (GAS5+GAS6+GAS7+GAS8 )/( 1 . 0+GAS9 ) 

100 MU=Z*1 .058E-06*16. 5273E00 

1000 FORMAT (/20X,48HWARNING! OUTSIDE OF VALIDITY RANGE OF CURVE FIT 
* r/ / 20X , 5HRHO = , 1PE15 . 8 , 5X , 3HT =,1PE15.8,/) 

RETURN 

END 
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ORIGINAL PAGE IS 
OF. POOR QUALITY 

SUBROUTINE UGAS2 (T, RHO, PR) 

INPUTS FOR SUBROUTINE : 

T = TEMPERATURE, IN KELVIN 
RHO = DENSITY, IN KG/M* *3 

OUTPUT : 

PR = PRANDTL NUMBER 

REAL T , RHO , PR , RHOO , Z , X , Y ,GAS1 ,GAS2 ,GAS3 , GAS4 
REAL GAS 5 ,GAS6 , GAS 7 ,GAS8 ,GaS9 
DATA RHOO/1 . 243/ 

X=T/1000 . 0E00 
Y=ALOG10 ( RHO/RHOO ) 

IF ( (Y.LT.-5.0E0O) .OR. ( Y .GT. 1 . 0E00 ) .OR . (X.GT. 1 . 5E01 ) ) 
* WRITE (6,1000) RHO , T 
IF (T.GT.500.E00) GO TO 10 
GAS1=7.16321E-01+1.1135E00*X 
GAS2=(5.58243E-06-7.16815E-05*X)*Y 
GAS3=(-7.72911E00+2.25827E-04*Y+1.44166E01*X)*X*X 
GAS4=(-1.47156E-07-2.28926E-07*X-2.88338E-08*Y)*Y*Y 
GAS5=-1 . 409 9E-01-3. 3505 5E-01*X 
GAS 6= (-2. 5597 5E -05+1. 5853E~04*X) *Y 
GAS7=(6.09194EOO-3.18345E-04*Y-1.32747E01*X)*X*X 
GAS8= ( 1 . 3742E-06-1.29479E-06*X+1.48302E-07*Y)*Y*Y 
GAS9=EXP(8.636-3.03E01*X) 

GO TO 100 

10 IF (T.GT.2.0E03) GO TO 20 

GAS1=6.766E-01+5.33391E-02*X 

GAS 2= (-2.01021E-02+4. 04905E~03*X) *X*X 

Z=GAS1+GAS2 

GO TO 110 

20 IF (T.GT.4.0E03) GO TO 30 

GAS1=5. 35204E-01+1.64262E-01*X 

GAS2=(-6.72637E-02+3.42314E-02*X)*Y 

GAS3=(-3.88497E-02-3.16248E-03*Y+3.05280E-03*X)*X*X 

GAS4=(-7.81832E-03+1.84389E-03*X-3.46855E-04*Y)*Y*Y 

Z=GA S1+GAS2+GAS3+GAS4 

GO TO 110 

30 IF (T.GT.6.5E03) GO TO 40 

GAS] =-2.39283EOO+1.28399EOO*X 
GAS: =(-7.675E-01+l .89502E-01*X)*Y 

GAS: =(-1.79581E-01-1.20286E-02*Y+8.30322E-03*X)*X*X 

GAS' =(-7.09301E-02 + 7.19471E-03*X-2.7837lE-03*Y)*Y*Y 

GAS: =3.06018EOO-1.20461EOO*X 

GkSt =(6.77677E-01-1.43868E-01*X)*Y 

GAS', = (1.62407E-01 + 7.85746E-03*Y-7.39086E-03*X)*X*X 

GAS8=(4.14157E-02-8.3635E-04*X-3.70369E-04*Y)*Y*Y 
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GAS9=EXP(-26 . 39E00+2 . 969EOO*X-5. 042EOO*Y-0 . 11 2E00*X*Y ) 
GO TO 100 

40 IF (T.GE.9.40E03) GO TO 50 

GAS1=6.13473E00-1 .54169EO0*X 

GAS 2= (1.08128E00-2. 04154E~01*X ) *Y 

GAS3=(1.43737E-01+9.91640E-03*Y-4.54467E-03*X)*X*X 

GAS4=(6.19987E-02-5.05808E-03*X+1. 56791E-03*Y)*Y*Y 

GAS5=~5. 4444 5E 00+1. 5845 9EOO*X 

GAS 6= (“1.10792E00+2.13203E-01*X)*Y 

GAS7= (-1 51000E-01-1. 00257E-02*Y+4 . 72964E-03*X) *X*X 
GAS8= (-7 80793E-02+7. 29918E-03*X-2.29357E-03*Y) *Y*Y 
GAS9=EXP ; 13 . 39E00-4 . 258EO0*X+2 . 298E00*Y-1 . 233E00*X*Y) 

GO TO 100 

50 IF (T.GT.ll . 5E03) GO TO 70 
IF (Y.GT.-2.5E00) GO TO 60 
GAS1=-3.24776E01+8.72772EOO*X 
GAS2=(-2.58872EO0+3.59O02E-01*X)*Y 
GAS3=(-7.61542E-01+2.10923E-02*Y+2.67953E-02*X)*X*X 
GAS4= (-1 . 9321E-01+8 . 94685E-02*X+5 . 64303E-02*Y ) *Y*Y 
GAS 5= 3 . 99935E01-9 . 68334EO0*X 
GAS6=(6.78337E00-9.07345E-01*X)*Y 
GAS7=(7.87932E-01-3.99108E-03*Y-2.64764E-02*X)*X*X 
GAS8=(6.97742E-01-1 .25709E-01*X-4.08833E-0?*Y)*Y*Y 
GAS9=EXP(105.8EOO-11.67EOO*X+31.67EOO*Y-3.: 3E00*X*Y) 

GO TO 100 

60 GASl=-2 .8075 5E01+6. 8040 6EOO*X 

GAS2=(-2.63243E00+4.0185E-01*X)*Y 
GAS3=(-5.45283E-01-1.63614E-02*Y+1.45424E-02*X) *X*X 
GAS4=(2.12026E-02-3.62386E-03*X+3.50018E-03*Y)*Y*Y 
GAS 5= 2 .82604E01-6.62279EO0*X 
GAS6=(2.06694E00-2.89135E-01*X)*Y 
GAS7= ( 5 . 26582E-01+1 . 13732E-02*Y~1 . 40944E~02*X ) *X*X 
GAS8=(-1 . 31445E-01+1. 50468E-02*X-9.13033E-03*Y)*Y*Y 
GAS9=EXP(-35.4lE00+2. 148EOO*X-1.481EOO*Y-0. 3704E00*X*Y) 

GO TO 100 

70 IF (Y.GT.-2.5E00) GO TO 90 

IF (T.GT.13.5E03) GO TO 80 
GAS1=6 . 08811E01-9. 88231 EO0*X 
GAS2=(9. 51872E00-9.95583E-01*X)*Y 
GAS 3= (5.48699E~01+2.67619E-02*Y-1. 03794E~02*X) *X*X 
GAS4=(5.0199E-01-2.55257E-02*X+8.45834E-03*Y)*Y*Y 
GAS5=-7.1667E01+1.2239E01*X 
GAS 6= (-1. 0995 9E 01+1. 2386 5E00*X)*Y 
GAS7=(-7.0869E-01-3 .61494E-02*Y+1.38445E-02*X)*X*X 
GAS8=(-3. 64368E-01+1 . 89712E-02*X+4 .04684E-03*Y) *Y*Y 
GAS9=EXP(-71 .89E00+2 . 248E00*X+0 . 4746E00*Y-0 . 9469E00*X*Y ) 
GO TO 100 

80 GAS1=2 . 9948 5E01-4. 1211 2EOO*X 

GAS2=(5.92879E00-5.27093E-01*X)*Y 
GAS3=(1.93397E-01+1 .19371E-02*Y-3.08939E-03*X)*X*X 
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GAS4= (4 .09472E-01-1. 78772E-02*X+9.49505E-03*Y)*Y*Y 
GAS 5= -2 .66557E01+3.05342EOO*X 
GAS6=(-S. 537 75E00 + 8.f 8359E-01*x: *Y 

GAS 7= (-9.53141E-02-1.98247E-02*''. +4 . 69853E-04*X) *X*X 
GAS8= (-7.622 32E- 01 +4.341 26E-02*X-5.10053E- 03* Y)*Y*Y 
GAS9=EXP(-540. 2E00+34. 3E00*X-14b.4EC0*Y+9.148E00*X*Y) 

GO TO 100 

90 GAS1 = ~3 . 1866 6E 00 + 8. 08818 E~01*X 

GAS2=(-4.00164E-01+3.59959E-02*X)*Y 

GAS 3= (-6.06519E-02-1 . 04205E-03*Y+1 . < 5243E~03*X) *X*X 
GAS4=(1.6658E-02-4.36487E-03*X-1.86593E-03*Y)*Y*Y 
GAS5=2.68501EOO-4.32123E-01*X 
GAS6= (1 .36103E-01+2. 5886E-02*X)*Y 

GAS7=(2.32842E-02-1.82391E-03*Y-4.09433E-04*X)*X*X 
GAS8=(-5. 37705E-02+1.0074E-02*X+2.05852E-03*Y)*Y*Y 
GAS9=EXP(-31.16E00+1.633E00*X+2. 395E00*Y-0.8707E00*X*Y) 

100 Z=GAS1+GAS2+GAS3+GAS4+ (GAS5+GAS6 K3AS7+GAS8 ) /( 1 . 0E00+GAS9 ) 

110 PR=Z 

1000 FORMAT( /20X, 48HWARNING! OU'SIDE OF VALIDITY RANGE OF CURVE FIT 
*,/,20X, 5HRHO = , 1PE1 5 . 8 , 5X, 3’IT =,1PE15.8,/) 

RETURN 

END 
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SUBROUTINE UGAS3 ( E , RHO , MU) 

INPUTS FOR SUBROUTINE : 

E = INTERNAL ENERGY, IN (M/S)** 2 
RHO * DENSITY, IN KG/(M**3) 

OUTPUT : 


MU » DYNAMIC VISCOSITY, IN KG/M-S 


5 


10 


20 


REAL E , RHO , MU , E0 , RHOO , Z , T , Y , F , GAS1 , GAS 2 , GAS3 , GAS 4 , 
*GAS5 , GAS 6 ,GAS7 , GAS8 , GAS9 
DATA RHOO, E0/1. 243, 78408. 4EOO/ 

Z=ALOG10 ( E/EO ) 

Y = ALOG 1 0 ( RHO/RHO 0 ) 

IF ( (Y.LT.-5.0E00) .OR. (Y.GT. 1 . OEOO) ) 

* WRITE (6,1000) RHO , E 
IF (Z.GT.0.44E00) GO TO 5 
T=0 . 4E00*E/287 . 06E00 

MU=1.462E-06*SQRT(T)/(1.0EOOai2.0EC3/T) 

RETURN 

IF (Z.GT.0.67E00) GO TO 10 

GASl=4.84547E-01+4 .67135E-01*Z 

GAS2=(5.71205E-04-1.43629E-03*Z)*Y 

GAS3=(2. 5511OE00-2 .33472E-04*Y-1.441O2EOO*Z)*Z*Z 

GAS4=(2 .534K;E-04-4. 72375E-04*Z+1 .86899E-05*Y)*Y*Y 

F=GAS1+GAS2+GAS3+GAS4 

GO TO 90 


IF (Z .GT. 1 . 75EOO ) GO TO 20 

GASl = -3. - : 1666E01 + 6.67883E01*Z 

GAS2=(-2. 43998E00+2. 12309E0O*Z)*Y 

GAS 3= (-3.69259E01-3. 08426E~01*Y+7 . 36486E00*Z ) *Z*Z 

GAS 4= ( -1 . 46446E-01+7 . 54423E-02*Z-2. 91464E~03*Y) *Y*Y 

GAS 5= 3 . 61757E01-6.11102E01*Z 

GAS 6= (2.40531EOO-2.05914EOO*Z)*Y 

GAS 7= (3.23911E01+2. 79149E-01*Y~5 .07640E00*Z)*Z*Z 

GAS8= (1 . 37916E-01-6. 72041E-02*Z42 .61987E-03*Y)*Y*Y 

GAS9=EXP( -3 .433E01-1 .823E00*Y+2 . 499E01*Z+6. E03E-01*Z*Y) 

GO TO 80 


IF (Z .GT. 2 . 5OEO0 ) GO TO 30 

GAS1=~1 .65147E02+2. 11028 E02*Z 

GAS 2= (-4.70948E00+2.78 >58E00*Z)*Y 

GAS 3= (~8.78308E01-1.28 i71E-01*Y+l . 27639E01* Z ) *Z*Z 

GAS4= (-3 . 19867E-01+1 . 7.H79E 01*Z + 3.86106E-03*Y) *Y*Y 

GAS5=2. 30407E02-2. 98055E02*/ 

GAS6=(-6.183O7E00+8.44595E0O*Z)*Y 

GAS7=(1.26933E02-2.61671E00*Y-1.' 7257E01*Z ) *Z*Z 

GAS6= (-2.30229E-02+2 . 25458E- 02* Z-4 . 41072E-03*Y) *Y*Y 

GASS=EXP( -6.882E01+8 .824EO0*Y+3.2O3EOl*Z-5.359E0O*Z*Y) 
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GO TO 80 

30 IF (Z.GT.2.85E00) GO TO 40 

GASl=-7 .09274E03+7 . 13648E03*Z 

GAS2=( -2 . 4601 4E0 2+1 . 65826E02*Z) *Y 

GAS3=(-2.37952E03-2.75487E01*Y+2.63465E02*Z)*Z*Z 

GAS4=(-3.49744EOO+1.28641EOO*Z-3.13711E-03*Y)*Y*Y 

GAS 5= 5 .26158E03-4.96701E03*Z 

GAS6= ( 2 . 03138E02-1 . 32984E02*Z) *Y 

GAS 7= (1.52424EC3 + 2. 15081E01*Y-1 . 504E 0E02*Z ) * 3*Z 

GAS8= ( 2 . 324 32EC 0-1 . 15997E00*Z+1 . 14862E-02*Y ) *Y*Y 

GAS9=EXP(-3 . 594E02-3 . 763E01*Y+1 . 319E02*Z+1 . 348E01*Z*Y) 

GO TO 80 

40 IF (Z.GT. 3 .15EG0) GO TO 50 

GAS1=-1.27748E03+1.29400E03*Z 

GAS2=(-3.60724E01+2.63194E01*Z)*Y 

GAS3= ( -4.22958E02-4 . 38228E00*Y+4 . 50571E01*Z) *Z*Z 

GAS4= ( -4 . 74425E-01+2 .89684E-01*Z+1.64048E-02*Y)*Y*Y 

F=GAS1+GAS2+GAS3+GAS4 

GO TO 90 

50 IF (Y.GT.-3.80E00) GO TO 70 

IF (Z.GT.3.19E00) GO TO 60 
GAS1=4.55919E03-4.21057E03*Z 
GAS2= (1. 03001 E0 1-2. 63478E01*Z)*Y 
GAS3= ( 1 . 29069E03+6 . 59587E00*Y-1 . 31413E02*Z) *Z*Z 
GAS4=(-8.28137E00+1 .9827E00*Z-1 .7287E-01*Y)*Y*Y 
F=GAS1+GAS2+GAS3+GAS4 
GO TO 90 

60 Z=E/E0 

GASl=-4 . 41792 E0 2+9. 7986E-02*Z 
GAS2=(-3.03148E02+7.6065E-03*Z)*Y 

GAS3= ( -5.5711E-05-3. 52836E-06*Y+8 .86148E-09*Z)*Z*Z 
GAS4 = ( -7.561E01-4. 76816E-04* Z-6 . 48859E00*Y ) *Y*Y 
GAS5=6.72387E04+3.28398E00*Z 
GAS6= ( 3 . 55009E04+2.72616EOO*Z)*Y 

GAS7=(2 . 13714E-03+3.42377E-04*Y-6.84897E-08*Z)*Z*Z 
GAS8= ( 6. 50886E03+3 . 8056E-01*Z+4 . 14116E02*Y ) *Y*Y 
GAS9=EXP(2 .978E01+5.415EOO*Y+1.713E-03*Z+3.115E-04*Y*Z) 

F=GAS1 +GAS 2+ GAS 3+GAS4 + ( GAS 5 +GAS6 +GAS7 +GAS8 ) / ( 1 . 0 -GAS 9 ) 

GO TO 90 

70 GASl=-6.4029E03+6.24254E03*Z 

GAS2= ( 1. 0327 9E02-8. 73181 E01*Z)*Y 
GAS3=(-2.02865E03+1.71878E01*Y+2.19907E02*Z)*Z*Z 
GAS4= ( -1. 22397 E0 1+3. 57830E00*Z-1. 27953E-01*Y ) *Y*Y 
F=GAS1+GAS2+GAS3+GAS4 
GO TO 90 

80 F=G£ S1+GAS2+GAS3+GAS4+ (GAS5+GAS6+GAS7+GAS8 ) /( 1 . 0+GAS9 ) 

90 MU=I . 74858 3E-05*F 

1000 FORMAT (/20X , 48HWARNING! OUTSIDE OF VALIDITY RANGE OF CURVE FIT 
* , /, 20X , 5HRHO = , 1PE15 . 8 , 5X , 3HE =,1PE15.8,/) 

RETURN 

END 
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SUBROUTINE UGAS4 ( E , RHO, K ) 

INPUTS FOR SUBROUTINE : 

E = INTERNAL ENERGY, IN (M/S)** 2 
RHO = DENSITY, IN KG/(M**3) 

OUTPUT : 

K = COEFFICIENT OF THERMAL CONDUCTIVITY, IN J/(KELVIN*M*S ) 

REAL E , RHO , K , EO , RHOO , T , Z , Y , F , GAS1 , GAS2 , GAS3 , GAS4 , GAS 5 , 
*GAS6 ,GAS7 ,GAS8 ,GAS9 
DATA RHOO ,E0/1 . 243EOO , 78408 . 4E00/ 

Z=AI OGIO (E/EO ) 

Y=AI OGIO ( RHO /RHOO ) 

IF ( (Y.LT.-5.0E00) .OR. ( Y.GT. 1 . OEOO ) ) 

*WRITE (6,1000) RHO, E 
IF ( Z.GT.0.44E00) GO TO 5 
T=0 . 4E00*E/287 . 06E00 

K=1 . 994E~03*SQRT (T)/ ( 1 . OEOO +11 2 . 0E00/T) 

RETURN 

5 IF (Z.GT.0.65E00) GO TO 10 

GAS1=1 .81003b9E-01+4.8126802EOO*Z 

GAS 2= ( -2 . 7231116E-02+1 .2691337E-01*Z)*Y 

GAS3=(-8 . 991 3 0 34EOO-1 .262408 5E- 01* Y+8. 96493 05E00*Z)*Z*Z 

GAS4= (-4 . 7198236E-03+9. 2328079E-03*Z-2 .9486327E-04*Y) *Y*Y 

F=GAS1+GAS2+GAS3+GAS4 

GO TO 200 

10 IF (Y.GT.-l.OOEOO) GO TO 130 
IF (Y.GT.-3.00E00) GO TO 70 
IF (Z.GT.1.25E00) GO TO 20 
GAS1=-1 . 0593 5E04+2. 3147 0E 04* Z 
GAS 2= (~7.41294E02+1.21724E03*Z)*Y 
GAS3=(-1 .67601E04-4.43184E02*Y+4.06631E03*Z)*Z*Z 
GAS4= ( 1 3510 5E 01 + 4. 94914E00*Z + 1 . 5538 5E00*Y)*Y*Y 
GAS5=1 .i)6032E04-2.31560E04*Z 
GAS 6= (7.46951E02~1.22465E03*Z)*Y 
GAS7=(1. 67604E04+4 . 45919E02 *Y-4 . 06258E03*Z ) *Z* Z 
GAS8=(-1 . 28615E01-5. 32398E00*Z-1 . 52956E00*Y) *Y*Y 
GAS9=EXP(-4. 219E01-4. 687E00*Y+2 . 812E01*Z+3 . 125EOO*Y*Z ) 
F=GAS1+GAS2+GAS3+GAS4+ (GAS5+GAS6+GAS7+GAS8 )/( 1 . 0-GAS9 ) 

GO TO 200 

20 IF (Z.GT.1.775E00) GO TO 30 
GAS1=3 . 7937 5E 03-7. 40351E03*Z 
GAS2= ( 3 . 29698E02-3 .55916E02*Z)*Y 
GAS3=(4.77122E03+1.00241E02*Y-1.00740E03*Z)*Z*Z 
GAS4= (1.97061E01-8.42554E00*Z+4.80494E~01*Y)*Y*Y 
GAS5=-4 . 5360 3E 03+9. 0560 5E03*Z 
GAS 6= (~4.95870E02+6.33563E02*Z)*Y 
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GAS7= ( -5. 95317E03-2 . 05442E02*Y+1 . 28945E03*Z ) *Z*Z 
GAS8=(-2.00087E01+1 .18851E01*Z-1.71735E-01*Y)*Y*Y 
GAS9=EXP(-3.318E01+3.158E-01*Y+1.863E01*Z-1.035EOO*Y*Z) 

GO r X> 190 

30 IF < Z.GT.1.93EOO) GO TO 40 

GAS1=2 .0665187 5E0 5- 3. 16564 5E05*Z 

GAS: =(-3.07322021E02+4. 57036377E02*Z)*Y 

GAS: =(1. 6182493 7E05-1 .555084 53E02*Y-2.7603957E04*Z)*Z*Z 

GAS 4 = (1. 9226026 5E00 _ 2. 2478809 4E00*Z~3. 0622601 5E~01*Y) *Y*Y 

GAS: =- 2. 06564 312E0 5+3. 18191 31 2E05*Z 

GASf =( 2.1 754228 5E03-2. 4 6670776E03*Z)*Y 

GAS' =( -1.63 597 062E05+7.16753174E02*Y+2.80926367E04*Z)*Z*Z 
GASi =( 3. 3952682 5E01 -7. 5384664 5E00*Z+1. 91 2143 71E00*Y)*Y*Y 
GAS' =EXP(-3.924E02-5.206E01*Y+2.054E02*Z+2.679E01*Y*Z) 

GO r iO 190 

40 IF ( Z.GT.2.60E0O) GO TO 50 

GAS1=7.1572625E04-9.2471625E04*Z 

GAS 2= (1.9646323E03-2.0280527E03*Z)*Y 

GAS? =( 3. 944610 5E04+4. 5673853E02*Y-5.5728672E03*Z)*Z*Z 

GAS4 = (-9.2131958E01 + l .2724541E01*Z-5.0568476E00*Y)*Y*Y 

GAS5=-3.2910781E04+4.2551211E04*Z 

GAS6=(1.4566331E03-2.2653745E03*Z)*Y 

GAS7=(-1.9476277E04+8 .4370288E02*Y+3.2389702E03*Z)*Z*Z 

GAS8=(-1 . 3324 594E02+1 .0591 53 3E02*Z+5. 86394 69E00*Y)*Y*Y 

GAS9=EXP(4.917E01+2.415E01*Y-2.455E01*Z-1.181E01*Y*Z) 

GO TO 190 

50 IF (Z.GT.2.69E00) GO TO 60 

GASI =1 .14 568 3E06- 1.23 752 5E06*Z 

GAS 2= (1.4024508E04-9.3467227E03*Z)*Y 

GAS 3= (4. 459305 6E 05+1. 533074E03M -5. 3608352 E04*Z)*Z*Z 

GAS 4= (2. 8485107E02-1. 096891 6E02* Z-l . 0955791E00*Y ) *Y*Y 

GAS5=-1.752087E06+1.79675E06*Z 

GAS6=(-1.3278737E05 + 9.8215562E04*Z)* { 

GAS7=(-6.0791744E05-1 . 811943E04*Y+6 . 7709875E04*Z )*Z*Z 

GAS8=(-1 . 3384084E03+5.2707324E02*Z+2.5904894E00*Y)*Y*Y 

GAS9=EXP(-1.798E02+7. 371E00*Y+6.73lEOl*Z-3.205EO0*Y*Z) 

GO TO 190 

60 GASl = -8 . 549962 5E 04 + 1. 173965 6E05* : 

GAS 2= (6.4563168E04-3.9551203E04* 0*Y 

GAS 3= (-4.8170254E04+6. 0816055E03 f Y+6 . 2052031E0: * Z ) *Z*Z 
GAS4=(2. 3473167E-01+1 .8871567E01 f Z+4. 0757723E00*Y) *Y*Y 
GAS5=5. 854688 3E04-9. 4 63487 5E04*Z 
GAS 6= (-6.6513812E04+4. 0899945E04-Z ) *Y 
GAS7=(4.2127227E04-6. 3717305E03*Y-5.7495195E03*Z)*Z*Z 
GAS8=(-1.0260344EOO-5.343277E01*Z-1.1017392E01*Y)*Y*Y 
GAS9=EXP(5.411E00+1.162E01*Y-1.082E00*Z-3.391E00*Y*Z) 
F=GAS1+GAS2+GAS3+GAS4+ (GAS5+GAS6+GAS7+GAS8 ) , ( 1 . 0-GAS9 ) 

GO TO 200 

70 IF (Z.GT.1.29E0O) GO TO 80 

GAS1=-1.22493E04+2.41071E04*Z 
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GASi!= (-1. 6182 9E 03 + 2. 2253 5E03*Z)*Y 
GAS3=(-1 .59261E04-7.53213E02*Y+3.53376E03*Z)*Z*Z 
GAS 4= (1. 9802 6E 00+5. 18483E00*Z+1.47851E00*Y)*Y*Y 
GAS5=1. 2248 6E 04-2. 41023E04*Z 
GAS 6= (1.61810E03-2.22571E03*Z)*Y 
GAS7=(1. 59235E04+7.53746E02*Y-3. 53168E03*Z)*Z*Z 
GAS8= ( -2 . 15482E00-5 . 05115E00*Z-1 ,48795E00*Y) *Y*Y 
GAS9=EXP(-3 . 111E01-4.444E00*Y+1. 944E01*Z+2.778E00*Y*Z ) 
F=GAS1+GAS2+GAS3+GAS4+ (GAS5+GAS6+GAS7+GAS8 )/( 1 . 0-GAS9 ) 
GO TO 200 

80 IF (Z.GT.1.85EOO) GO TO 90 

GAS1=3.18060E03-6.69664E03*Z 

GAS 2= (4.33382E01-2.14649E02*Z)*Y 

GAS 3= (4. 41377E03+9.41359E01*Y~9. 29758E02*Z ) *Z*Z 

GAS 4= (-3.62190E01+1 . 15538E01*Z-2 . 14621EOO*Y) *Y*Y 

GAS5=-5.98764E03+1 .29243E04*Z 

GAS 6= (~2.72261E02+5.42378E02*Z)*Y 

GAS7=(-9.03293E03-2.11787E02*Y+2.07831E03*Z)*Z*Z 

GAS8=(2.74179E01-5.68578EOO*Z+1.91217EOO*Y)*Y*Y 

GAS9=EXP(-1 .854E01+7.11EOO*Y+1.068EOl*Z-5.449E00*Y*Z) 

GO TO 190 

90 IF ( Z.GT.2.OEO0) GO TO 100 

GAS! =5 . 14024E04-7 . 52733E04*Z 

GAS 2= (-3.30889E02+3.11550E02*Z)*Y 

GAS 3= (3. 6653 9E 04-7. 41227E01*Y~5. 9301 5E03*Z)*Z*Z 

GAS4=(-4.84164E01+2.23133E01*Z-9.19118E-01*Y)*Y*Y 

GAS5=-1.80898E05+2.82532E05*Z 

GAS 6= (-1.01053E03+9.75576E02*Z)*Y 

GAS7=(-1.47220E05-2.33631E02*Y+2.55940E04*Z)*Z*Z 

GASf=(3.2868lE0O-1.76588E00*Z-1.54962E-Ol*Y)*Y*Y 

GASS =EXP(-4 . 104E01 + 6. 507E01*Y+2 . 083E01*Z-3 .472E01*Z*Y) 

GO TO 190 

100 IF { Z.GT.2.58E00) GO TO 110 

GAS1=5 .1131824E04- 6 . 664875E04*Z 

GAS 2= (2.02171E03-J . 9306292E03*Z ) *Y 

GAS 3= (2.8762395E04+4.3353467E02*Y-4.1064609E03*Z)*Z*Z 

GAS 4= (-&.4970047E01+1 . 7925919E01*Z-6 . 2576542F 00*Y ) *Y*Y 

GAS5=-6. 276815 6E 04+8. 601587 5E04*Z 

GAS 6= ( -1. 00020 36E0 3+6. 2 537280 E02*Z) *7 

GAS7= (-3. 95782 7E04-3. 8467 377E01* Y+6. 12953E03*Z)*Z*Z 

GAS8= (-1.0591702E02+7. 63614 2E01* Z+5. 938859E00*Y) *Y*Y 

GAS9=EXP(-3.901E00+2.418E01*Y+1. 374E00*Z-1.145E01*Y*Z) 

GO TO 190 

110 IF (Z.GT.2.73E00) GO TO 120 

GAS1=1. 0088046E06-1. 08632 1E06*Z 

GAS2=(1. 3844801E04-9.7268516E03*Z)*Y 

GAS 3= ( 3. 898532 5E05+1. 709166 5E 03* Y- 4. 6621066E04*Z ) *Z*Z 

GAS 4= (1. 484072 6E02-5.2645004E01*Z-1. 5477133E~01*Y) *Y*Y 

GAS5=-1.073351E06+1.14571E06*Z 

GAS 6= (-1.9343957E04+1.3366211E04*Z)*Y 
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GAS 7= (-4. 0670987 E05-2. 2955198 E03*Y+4.7999871E04*Z)*Z*Z 
GAS8=( -4.1 01 6724E02+1. 49941 48E02*Z-1.9779787E00*Y)*Y*Y 
GAS9=EXP(-1 . 026E02+6. 302E01*Y+3 . 819E01*Z~2 . 431E01*Y*Z ) 
GO TO 190 

120 GAS1 = -9.6638500E04 + 1.3206488E04* Z 

GAS2=( -4. 74 5810 5E04+2. 359687 5E04*Z)*Y 

GAS3= (1 .8602773E04-2 . 306802E03*Y -4 . 0413552E03*Z ) *Z*Z 

GAS4= (- 5. 3564258E03+2 . 2433904E03*Z+2 . 5188145E02*Y) *Y*Y 

GAS 5=1 .0962581E05-2. 99011 6E04*Z 

GAS6= ( 4 . 7883496E04-2 ,3785383E04*Z)*Y 

GAS 7= (“1. 175396 9E 04+2. 290552 2E03*Y+3. 130439 9E03*Z)*Z*Z 

GAS 8= (5.473418E03-2.3208018E03*Z -2 . 6570068E02*Y ) *Y*Y 

GAS9=EXP(-3.107E01+1.082E01*Y+1.047E01*Z-3.047EOO*Y*Z) 

F=GAS1+GAS2+GAS3+GAS4+ (GAS5+GAS6+GAS7+GAS8 )/( 1 . 0-GAS9 ) 

GO TO 200 

130 IF (Z.GT.1.40E00) GO TO 140 

GAS1=-1 . 58386E03+3 ,49223E03*Z 

GAS 2= (~8.39834E02+1.09565E03*Z)*Y 

GAS3=(-2.56175E03-3.56197E02*Y+6.25145E02*Z)*Z*Z 

GAS4=(-1.22407E01+7.65634E00*Z+2.58235E-01*Y)*Y*Y 

GAS 5=1 . 5802 5E 03~3. 47664E03*Z 

GAS 6= (8.39588E02~1.09490E03*Z)*Y 

GAS7=(2.54682E03+3.55674E02*Y-6.18504E02*Z)*Z*Z 

GAS8= (1 .20843E01-7 . 44857E00*Z -2 . 91202E-01*Y ) *Y*Y 

GAS9=EXP(-2.171E01-4.342EOO*Y*-1. 316E01*Z + 2.632EOO*Y*Z) 

F=GAS1+GAS2+GAS3+GAS4+ (GAS5+GAS6+GAS7+GAS8 )/( 1 . 0-GAS9 ) 

GO TO 200 

140 IF (Z.GT.1.91E00) GO TO 150 
GAS1=7.89255E02-1 .91743E03*Z 
GAS 2= (3.59227E02-4.44070E02*Z)*Y 
GAS3=(1.39463E03+1.34083E02*Y-3.13446E02*Z)*Z*Z 
GAS4= (1 .90681E01-1 . 09285E01*Z+4 . 24933E~02*Y) *Y*Y 
GAS5=-1.31401E03+3.13134E03*Z 
GAS6=(-5.18755E02+6.80268E02*Z)*Y 
GAS7=(-2.32493E03-2.21393E02*Y+5.52563E02*Z)*Z*Z 
GAS8=(-3.32001E01+2.11819E01*Z-4.75163E-01*Y)*Y*Y 
GAS9=EXP(-5.025EOl-8.412E0O*Y+2.982EOl*Z+3.509E00*Y*Z) 
GO TO 190 

150 IF (Z.GT.2.05E00) GO TO 160 
GAS1=3.58691E04-5.16852E04*Z 
GAS2=(-6.30189E02+6.63314E02*Z)*Y 
GAS 3= (2.47471E04-1 . 73538E02*Y~3 . 93167E03*Z ) * Z*Z 
GAS4=(-4.23871E01+2.08048E01*Z-1.05512EOO*Y)*Y*Y 
GAS5=-1.10522E05+1.67591E05*Z 
GAS6=(4.61877E03-4.94930E03*Z)*Y 
GAS ’=(-8.46558E04 + 1.32441E03*Y+1.42438E04*Z)*Z*Z 
GASf =(2.25065E01-1.10316E01*Z+9.62887E-01*Y)*Y*Y 
GASS =EXP(-1 .68lE02 + 7.063E01*Y+8.75E01*Z-3. 75E01*Y*Z) 

GO TO 190 

160 IF (Z.GT.2.57E00) GO TO 170 
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GAS1 =3.1899562E04~4.2186664E04*Z 

GAS2 =( 2. 305560 3E0 3-1 . 98970 17E0 3* Z)*Y 

GAS 3= (1.849998E04+4. 256181 6E02*Y-2. 6808696E03*Z ) *Z*Z 

GAS4=(-1.6195114E01+5.8640623EOO*Z-3.6172504EOO*Y)*Y*Y 

GAS 5= ~5 ,7594039E04+7.9328437E04*Z 

GAS 6= (-1. 927598 9E 03+1. 673054 4E03*Z)*Y 

GAS7 =(-3.6473008E04-3.6100732E02*Y+5. 597543E03* Z ) *Z*Z 

GAS8=(-7.920808E01+4.0542084E01*Z+2.1495867EOO*Y)*Y*Y 

GAS9=EXP(-5.733E01+2. 088E01*Y+2. 592E01*Z-9.793E00*Y*Z) 

GO TO 190 

170 IF (Z.GT.2.75E00) GO TO 180 

GAS1=7 .0838087E05-7. 561991 9E05*Z 

GAS 2= ( 3. 9503091E03-2. 7381802 E03*Z)*Y 

GAS3=(2 .68881 81E0 5+4 . 7728687E 02*Y-3 . 183816E04*Z )*Z*Z 

GAS4= (-1.2532251E02+4. 773478 7E01*Z- 4 . 0148029E00*Y) *Y*Y 

GAS 5= -2 . 521632 5E05+2. 172776 9E05*Z 

GAS6=(9.2882383E03-7.780918E03*Z)*Y 

GAS 7= ( ~5. 6539297E0 1+1. 6120212E03*Y+3. 941924 8E03*Z)*Z*Z 

GAS8=(1.8537296E02-7.1010757E01*Z+1.1307096E00*Y)*Y*Y 

GAS9=EXP(-1.786E02+2.18E-01*Y+6.714E01*Z-4.739E-01*Y*Z) 

GO TO 190 

180 GAS! =3. 185503 7E0 5-3. 30411 56E05*Z 

GAS2=(2.2983352E04-1.6623461E04*Z)*Y 
GAS:=(1.13848E05+3.0098223E03*Y-1.3020133E04*Z)*Z*Z 
GAS4 =( -1.8 5990 39E02+6. 984068 3E01*Z-7.7371645E00*Y)*Y*Y 
F=GAS1+GAS2+GAS3+GAS4 
GO TO 200 

190 F=GAS1+GAS2+GAS3+GAS4+ (GAS5+GAS6+GAS7+GAS8 )/( 1 . 0+GAS9 ) 

200 K=1 . 8791 5E-02*F 

1000 FORMAT(/20X, 48HWARNING! OUTSIDE OF VALIDITY RANGE OF CURVE FIT 
* ,/ , 20X, 5HRHO = , 1PE15 . 8 , 5X, 3HE =,1PE15.8,/) 

RETURN 

END 
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